#ifndef LATTICE_HPP
#define LATTICE_HPP
#include "Types.hpp"
#include "Grid.hpp"
#include <cmath>
#include <map>
#include "VTKFileWriter.hh"

class Lattice
{
    public:
        Lattice(uint lengthx,uint lengthy,real &omega);
		Lattice(real &omega, std::string &pgm_file);
		Lattice();
        ~Lattice(){}
        inline real fequilibrium(uint x,uint y, uint f);
        void initialize(real density, real velX, real velY);
        inline const uint& getLengthx() {return lengthx_;}    
        inline const uint& getLengthy() {return lengthy_;}
        inline const real& getOmega() {return omega_;}
        inline void setOmega(real nOmega){omega_=nOmega;}
        inline void updateVelocity(uint x, uint y);
        inline void updateDensity(uint x, uint y);
        inline bool isFluid(uint x, uint y){return (nodesType_(x,y)==FLUID)?true:false;}
        void streaming(real vwx,real vwy);
        void collision();
        void streamingCollision(real vwx,real vwy);
        inline void swapPDFsGrids(){/*swap(src_,dest_);*/}
        void printF(const std::string &vtkFile);
        void copyGrid();
        void comparePDFsGrids();
        void checkMaps();

    private:
        uint lengthx_;
        uint lengthy_;
        real omega_;
        PDF_Field dest_;
        PDF_Field src_;
        V_Field vs_;
        D_Field rhos_;
        Flags nodesType_;
        std::map<uint,uint> inversePDF;
        std::map<uint,real> weights;
};

inline real Lattice::fequilibrium(uint x,uint y, uint f)
{
    real dotProduct=vs_(x,y,0)*eX[f]+vs_(x,y,1)*eY[f];

    real sum = weights[f]*rhos_(x,y)*(1.0+3.0*dotProduct+4.5*dotProduct*dotProduct-1.5*(vs_(x,y,0)*vs_(x,y,0)+vs_(x,y,1)*vs_(x,y,1)));

    return sum;
    //return weights[f]*(rhos_(x,y)+3.0*dotProduct+4.5*dotProduct*dotProduct-1.5*(vs_(x,y,0)*vs_(x,y,0)+vs_(x,y,1)*vs_(x,y,1)));
}

inline void Lattice::updateDensity(uint x, uint y)
{
    real sum=0;
    for(uint k=0;k<DIMEN;++k)
        sum+=dest_(x,y,k);
    rhos_(x,y)=sum;
}

inline void Lattice::updateVelocity(uint x, uint y)
{
    real sumX=0, sumY=0;
    for(uint k=0;k<DIMEN;++k)
    {
        sumX+=dest_(x,y,k)*eX[k];
        sumY+=dest_(x,y,k)*eY[k];
    }
    vs_(x,y,0)=sumX/rhos_(x,y);
    vs_(x,y,1)=sumY/rhos_(x,y);
// if (vs_(x,y,0) >10 ||vs_(x,y,1) >10 )
// std::cout << vs_(x,y,0) << " " << vs_(x,y,1) << std::endl;

}

#endif //LATTICE_HPP
